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Oscillatory motion of sheared nanorods beyond the nematic phase 
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We study the role of the control parameter triggering nematic order (temperature or concentra- 
tion) on the dynamical behavior of a system of nanorods under shear. Our study is based on a set of 
mesoscopic equations of motion for the components of the tensorial orientational order parameter. 
We investigating these equations via a systematic bifurcation analysis based on a numerical contin- 
uation technique, focusing on spatially homogeneous states. Exploring a wide range of parameters 
we find, unexpectedly, that states with oscillatory motion can exist even under conditions where 
/--v ' the equilibrium system is isotropic. These oscillatory states are characterized by wagging motion of 

^sj , the paranematic director, and they occur if the tumbling parameter is sufficiently small. We also 

present full non-equilibrium phase diagrams, in the plane spanned by the concentration and the 
t^*V shear rate. 
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I. INTRODUCTION 

The non-equilibrium dynamics of fluids composed of interacting, rod-like particles under external shear flow is a 
long-standing problem receiving continuous attention both from a fundamental perspective and in applications [lHj] . 
Renewed theoretical interest results from the discovery of stable, oscillatory orientational modes as well as rheochaos 
under conditions where the particles strongly interact [5], Q . This is typically the case at high densities and/or low 
temperatures, depending on whether the system is rather lyotropic (density-driven) or thermotropic (temperature- 
driven). In thermal equilibrium, strong coupling between theparticles leads to various mesophascs with long-range 
orientational ordering and different degrees of positional order p], [8| . Under shear, one observes not only shifts of these 
various transitions [9(, but also oscillatory motion (and other dynamical modes) of the entire nematic director. Such 
behavior is predicted, e.g. by various continuum theories (such as those of Doi and Hess) [KH13I for the dynamics of 
the orientational order parameter of the system, as well as in particle-based computer simulated techniques such as 
(event-driven) Brownian dynamics and multi particle collision dynamics |14| - [l6j . Experimentally, oscillatory states 
have been seen, e.g., in suspensions of fd- viruses in plane Couette flow geometries |17lll8|. Again, the corresponding 
equilibrium system (at zero shear rate) is a nematic. 

However, despite the large number of studies already carried out, there are still some rather fundamental questions 
left. One of these is whether oscillatory motion can occur only if the system in the stable or metastable nematic 
regime, as it is often assumed. In other words, is it possible that spontaneous oscillatory motion develops on the 
isotropic side of the shear-distorted phase diagram? And if such oscillations indeed exist, how do they transform into 
other dynamical states if the shear rate or other system parameters are changed? Or, in the language of non-linear 
dynamics, what is the nature of the underlying bifurcations? 

In the present paper we aim to shed light on these issues on the basis of the mesocopic approach suggested by Hess 
[10l - ll2[ | and later used in many subsequent studies J19l422| . However, contrary to these earlier studies we here employ 
a systematic bifurcation analysis, which allows us to explore a much wider range of parameters and, in addition, to 
fully characterize the nature of the dynamical transitions. 

In the Hcss-Doi approach, the relevant dynamic variable is the tensorial orientational order parameter whose 
dynamics is determined by the (Couette) shear, on the one hand, and relaxational contributions (derived from a 
Landau free energy), on the other hand. By using a suitable tensor basis, the corresponding equation of motion 
transforms into a set of five, non-linear first-order differential equations for the expansion coefficients. It is well 
established that these equations predict rich state diagrams involving a variety of oscillatory states and even chaos. 
Here we analyze the dynamical behavior on basis of a numerical path continuation method, exploring two bifurcation 
parameters, that is, the shear rate and a coupling constant). Previous works by Alonso et al. [23| and Forest et al. 
[24j also used bifurcation analysis. However, they limited the description to orientations within the shear plane [231 ] or 
focused on the behavior for weak shear j24|. Another bifurcation analysis was presented by Rey [25|], who considered 
(contrary to our work) a system under steady biaxial stretching flow. 

Here we consider a steady uniaxial shear flow, but do not restrict the resulting orientational dynamics. Our bifurcation 
analysis leads to results that indicate there is indeed an "island" in parameter space where the equilibrium system 
is a stable isotropic state, whereas shear induces an oscillatory (wagging) motion. Moreover, we can extract entire 
non-equilibrium state diagrams revealing a close relationship to experimental results. 

The paper is organized as follows. We first we give an overview of the model system and the governing dynamical 
equations. Next we apply the bifurcation analysis to a specific state parameter within the nematic phase, the main 
purpose being to compare with earlier results obtained by simple numerical integration. We then move towards 
previously unexplored parameter regions and set up full non-equilibrium diagrams. Finally, we give a conclusion and 
an outlook. The paper is rounded up by two appendices describing details of our method(s). 

II. THEORY 

A. Second rank alignment tensor 

We are interested in the dynamics of a suspension of rodlike particles under shear. As a method of investigation we 
employ a mesoscopic approach originally suggested by Doi and Hess PJj-QjJ [26[ where the basic dynamical variable is 
the orientational order parameter a(r,i) (in the literature also often Q). The latter describes the average orientation 
of the rods in a small, yet macroscopic volume of space. The alignment tensor is linked to the electric susceptibility 
tensor [27[, and can be measured optically by birefringence. In an unsheared (i.e., equilibrium) nematic system it is 
usually sufficient to describe the orientational order by a vector, the so-called director, which indicates the average 
direction of the particles. The strength of this average orientation is then described by the Maier-Saupe parameter S. 
However, in the presence of shear the orientational order can become biaxial |19l.l28j. Therefore, the order parameter 



a(r,t) must be a second-rank tensor, which we define according to 
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a(r,t) = J — (u<g> u- -Tr(u <g> u)) ox = \l — { u®u ) or , (2.1) 

where the unit vector u = u(r, t) stands for the orientation of the rod at position r at time t, ® denotes the dyadic 
product, and Tr(. . .) denotes the trace of a matrix. As it follows from the first member of Eq. (|2.ip . the tensor a(r, t) 
is symmetric and traceless (which we henceforth indicate by using the symbol ttt ). Moreover, the average (. . .) or 
appearing in Eq. (|2.1|) is defined as 

(...)«= / d 2 u...p OI (u,r,t), (2-2) 
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involving the orientational distribution function p° r {u, r, t) [27|. The integral in Eq. (J2.2I) is performed with respect to 

the angles describing the vector u. The orientational distribution is defined as p° T (u, r, t) = AT - 1 Q^j = i 3( u ~ u i(t)))ens, 
where Ui is the microscopic orientation of particle i (i — 1, . . . , N), and (. . .) ens is an ensemble average in a small 
volume dV around the space point r at time t. Combining Eq. (|2.2|l with the last member of Eq. (|2.1|) . we see that 
the order parameter tensor a(r, t) corresponds to the second moment of p or (u, r, t). In other words, the orientational 
distribution function provides the link between the mesoscopic quantity a(r,t) and the microscopic degrees of freedom, 

Ui. 

The eigenvalues of the tensor a are pi,p2,P3 and they are related to the eigenvectors ("principal axes") l,m,n. The 
tensor a can then be written in the form a = p,\l ®l + pirn <S> rn + p^n ® n [23. The principal axis related to the 
largest eigenvalue is considered as the director of the system. Commonly, one assumes that this role is played by ps, 
such that the director is given by n. For uniaxial order, n is the only relevant axis, and p% is linked to the Maier-Saupe 
parameter S via ^3 = 2^/5/65 [301, where S = (Pzfa ■ n)) or and Pi denotes the second Legendre polynomial. For 
the mo re g eneral case of biaxial order, the principal axes can be conveniently visualized using super-quadratic tensor 

glyphs mm. 

B. Equations of motion for the homogeneous sheared system 

In the absence of shear, the ordering behavior of the system is essentially controlled by either the concentration, 
c, or the temperature, T. Specifically, c is the control variable for lyotropic liquid crystals, and for most colloidal 
rods (e.g. fd- viruses [131), whereas thermotropic liquid crystals are controlled by the temperature. Nematic phases 
typically occur at high concentrations or low temperatures, respectively. In the present paper we describe lyotropic 
and thermotropic systems on the same footing in terms of an effective, dimensionlcss "temperature" 9 defined as 

l-T*/T 



1 - T*/T K 
1 - c/c* 



(thermotropic), 
(lyotropic). (2.3) 



1 - c K /c* 

In Eqs. (12.31). 7 V and ck correspond to the temperature and concentration where the isotropic and nematic phases 
coexist [271 . |29| . In other words, Tk and ck define the location of the first-order isotropic- nematic transition of the 
system. These parameters are sometimes also referred to as " clearing point" due to the related change of the optical 
properties of the system. From Eqs. (|2.3[) it follows that at T = Tk (c = ck), the effective temperature 6 — 1. On the 
other hand 6 — at the " pseudocritical" temperature T* or pseudocritical concentration c* . For T < T* (c > c*) the 
isotropic phase is globally unstable. Between the lower spinodal temperature and the clearing point temperature Tk 
(clearing point concentration Ck), i-C < 6 < 1, without flow the nematic phase is stable and the isotropic phase 
is metastable. In the range between temperature Tk (concentration ck) and T u (c u ), i.e. 1 < 6 < 9/8, the nematic 
phase in the absence of shear is metastable and the isotropic phase is stable. That is, T u (c u ) marks the upper limit 
(lower limit) of the metastable nematic phase [33J|. The stability of the two phases and thus, the equilibrium behavior 
of the system, is determined by the (dimensionless) Landau-de Gennes (LG) free energy [2[ 

1 

$(a) = -Tr(a-a) - \/6Tr(a-a-a) + -(Tr(a-a)) 2 , (2.4) 

where we have focused on spatially homogeneous states (i.e., no gradient terms). Also, we have defined the rcscalcd 
order parameter 

a=— , (2.5) 
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where ax is the value of the parameter a defined as 



y/Tr(a ■ a). (2.6) 



evaluated at coexistence, that is, at 7k or ck- We note that, for the special case of uniaxial orientational order, the 
parameter a is proportional to the Maier-Saupe order parameter, that is, a = ^/bS. 

In non-equilibrium the LG free energy, or rather its derivative with respect to the order parameter, governs the 
relaxational dynamics of the system. This derivative follows from Eq. (|2.4[) as 

$'(a) = — = 8a- 3^6 a~a + 2Tr(aa)-a. (2.7) 

da 

In the presence of shear flow, the relaxational dynamics competes with the dynamics induced by flow field v(r). The 
resultin g te rms entering the equation of motion for a(t) can be derived from a generalized Fokker-Planck equation 
[111 . Il3l . |26| or, alternatively, from irreversible thermodynamics [lOj. In these equations, the influence of shear is 
captured by the symmetric and antisymmetric part of the velocity gradient 

r^((Vu) T + Vu) (2.8) 

n=i((Vt>) T -VtO, (2.9) 

where T and O are the so-called deformation and vorticity, respectively. Here we focus on the case of plane Couette 
flow characterized by a linear velocity profile v = / jye x , where 7 is the shear rate and e x is a unit vector. Thus, the 
shear plane is spanned by the x- and y-directions and the unit vector e z is orthogonal to the shear plane. 

Combining relaxational and shear-induced terms, the equation of motion for a homogeneous system (i.e., a(r,t) = 
a(t)) in reduced units reads |19j 

^ = 2fTa +2aV^a - & {a) + J ^ \ K T . (2.10) 



In Eq. (|2.10[) , the shear flow enters through the three terms involving the quantities T and ft defined in Eqs.(| 
respectively. Specifically, the first term of the right side of Eq. (|2.10[) describes the impact of the flow vorticity fi, 
while the second term couples the deformation rate T linearly to the alignment tensor. The third term represents 
the relaxational part determined by the LdG free energy [see Eq. (|2.7[) ]. Finally, the last term on the right side of 
Eq. (|2.10[) involves the so-called tumbling parameter Ak, which determines the impact of the external perturbation 
due to the flow. This coupling parameter can be related to the axis ratio q = L/D (with L length, D width) of the 
particles via [llj 

V5a K <r + 1 

Specifically, one has Ak = for spherical particles whereas Ak > for elongated particles. As an example we briefly 
consider fd- virus particles which have been experimentally studied under shear in [17|, LL8( . These particles have a 
rather large aspect ratio of q = L/D « 130 [3j,|35|, such that the ratio (q 2 — l)/(q 2 + 1) ~ 1. To determine the value 
of ok, we use the fact that typical values of the Maier-Saupe parameter at coexistence are about S w 0.6 — 0.7 [34j | 
and a = y/^S. Inserting these values into Eq. (|2.11[) we obtain A^ w 0.6 — 0.7. 

It is common to transform the tensorial equation (|2.10[) into a system of scalar equations 13611 . This is done by 
expanding a and the other tensors appearing in Eq. ([2.10J1 into a tensorial basis set (see, e.g., |37j). yielding 

«o = - 0o - -%/3o"Ya 2 



di = — (j)\ + ja 2 
a 2 = - 



- 7°i + -VsXrJ - -^30-700 



d3 = - 03 + -H* 7 + !) a 4 

a 4 = - 04 + tt7(^ - 1)03, (2.12) 



where 

00 =(0 - 3a + 2a 2 )a + 3(a 2 + a 2 ) - -(a 2 + a 2 ) 

0i =(6> + 6a + 2a 2 )ai - -\/3(a^ - a\) 

02 =(9 + 6ao + 2a 2 )a2 — 3v3a3a4 

03 = (9 — 3a + 2a 2 )a 3 — 3v3(aia3 + 0204) 

04 = (9 — 3ao + 2a 2 )a4 — 3v3(a2Q,3 — 0104). (2-13) 

In Eqs. (|2.13p . a 2 = Yli=o a f • I n Appendix 1X1 we present some general remarks concerning the structure of Eqs. (|2 . 12[) 
and (|2.13[) and the consequences for the dynamic behavior of the system. 

In previous studies [381 ] it has been shown, for the case of planar Couette flow, the parameter a has only minor 
importance. Therefore, we henceforth set a = 0. 

The coupled set of ordinary differential equations (ODEs) given in Eqs. (|2. 12[) can be studied using methods from 
nonlinear dynamics. Here we perform a bifurcation analysis using numerical continuation techniques, specifically, the 
software package MATCONT j39| . Some details of this method can be found in Appendix [B] 

III. NUMERICAL RESULTS 

The subsequent section is divided into three parts. In the first part (Sec. IIII A[) . we aim to validate the results of 
our bifurcation analysis by comparing with well-established results from direct numerical integration [19| | . To this end 
we consider the temperature 9 = and calculate a state diagram in the plane spanned by Ak and 7. The temperature 
9 = corresponds to the lower spinodal temperature T* (upper spinodal concentration c*). In other words, for 9 < 
without shear the nematic phase is the only stable phase, while the isotropic phase is (locally and globally) unstable. 
For 9 = a large body of results for this specific value already exists [ijj, [20, [23, [30] ■ Beyond the pure numerical com- 
parison, however, we also analyze in detail the nature of the observed bifurcations which are, so far, largely unknown. 
In the second part (Sec. IIIIBp we apply the bifurcation analysis to systems at larger values of 9 corresponding to the 
isotropic side of the equilibrium I-N transition [see text below Eqs. (|2.3[) ] . In particular, we demonstrate the existence 
of oscillatory solutions. As a summary and overview of our results, we finally present in Sec. IIII CI full non-equilibrium 
phase diagrams (in the plane 9-j or 07, respectively) for selected values of the tumbling parameter. Later in this 
section we also provide codimension-1 diagrams that help to interpret the codimension-2 state diagrams. 

A. Nematic Phase 

In this subsection we set 9 = 0, where the unsheared equilibrium system is in the nematic phase. Switching on the 
shear flow the system can be brought into a variety of time-dependent and steady states characterized by a particular 
behavior of the order parameter a(t) [19|, |28[. Here we investigate the dynamics in the plane spanned by the tumbling 
parameter Ak and the shear rate 7. 

In Fig. [1] we show a composite of data obtained by direct numerical integration (colored areas), on the one hand, 
and the bifurcation analysis (curves), on the other hand. Within the calculations based on direct integration, we 
have considered a narrow grid of parameter pairs (Ak , 7) ■ The resulting dynamics of the alignment tensor was then 
automatically classified (ignoring transient structures) and colored accordingly. The curves in Fig. [T] were obtained 
using the codimension-2 bifurcation analysis described in the Appendix. 

The different regions appearing in the phase diagram at 9 = are the wagging state (W), the tumbling state (T), 
kayaking- wagging (KW), kayaking-tumbling (KT), and the (flow-) alignment state (A). Apart from the last (stationary) 
state, where the director is "frozen" along a direction within the shear plane, all other states mentioned so far are 
characterized by a time-dependent, oscillatory behavior of the coefficients a,(t) of a(t). Physically, these oscillations 
correspond to oscillations of the nematic director either within the shear plane (W, T), or out o/the shear plane (KW, 
KT). In addition to these regular states, the diagram at 9 = also involves complex or chaotic (C) states, as has been 
shown in previous investigations [19|, |28|. In the present study, however, we will rather focus on the regular states. 

We note that the colored areas in Fi g, p] which have been obtained by direct integration, agree quantitatively 
with earlier numerical results (see, e.g., Qli). More importantly, however, we can see from Fig. [T] that these colored 
areas are closely bounded by the curves stemming from our bifurcation analysis. This already indicates that the two 
approaches yield consistent results for the general dynamical behavior of the system. In the following wc discuss in 
more detail the boundaries between the states, that is, the nature of the underlying bifurcations. 




FIG. 1. (Color online) Dynamical state diagram in the plane spanned by the tumbling parameter Ak and the shear rate 7 at the 
reduced temperature 9 = (and a = 0). Areas are computed via direct numerical integration of the ODEs [see Eqs. (|2,12|l and 
(12.121) ]. whereas the curves are found from a codimension-2-bifurcation analysis. The areas are colored and labeled according 
to the dynamical state (A= Alignment, W= Wagging, T=Tumbling, KW=Kayaking Wagging, KT=Kayaking Tumbling, and 
C=Complex/Chaotic). The bifurcation branches are labeled as LP=Limit Point (saddle-node), H=Hopf, PD=Period Doubling, 
LPC=Limit Point of Cycles, BPC= Branch Point of Cycles, and NS = Neimarck Sacker. In addition, the diagram involves 
some special codimension-2 bifurcation points. These are the CP= Cusp Point, CPC= Cusp Point of Cycles, ZH= Zero Hopf, 
and BT= Bogdanov Takens point. 



From the steady state to in-plane oscillations 



We begin by considering the situation encountered at high shear rates and large tumbling parameters (see upper 
right corner of the diagram Fig.Q]). In this case, the system settles into a "flow- aligned" (A) state where the nematic 
director remains temporarily fixed along a direction within the shear plane. From a physical point of view, this A state 
may be further characterized by the magnitude of the nematic ordering under shear and the flow angle tp between 
director and flow direction. Mathematically, the A state corresponds to a (stable) fixed point of the dynamical system 
Eq. (|2.12[) . In the case of 6 = considered in this Section, there is only one stable fixed point and thus, only one 
type of alignment. More generally, however, the A state can also involve more than one fixed point. An example 
will be given in Sec. IIIIBI where we consider a higher temperature 9 and find, within the A state, a coexistence of 
paranematic and nematic flow alignment [4fJ. 

Starting from the A state, we now decrease the tumbling parameter Ak- Provided that the shear rate is sufficiently 
large (7 > 4.5), the system then encounters a dynamical transition into the wagging (W) state (see Fig. [T]). In the W 



state the angle ip between the director (which still lies in the shear plane) and the flow direction oscillates periodically 
between a minimal and a maximal value. This angular motion is accompanied by an oscillation of the magnitude of 
nematic order. 

The dark green line (labeled by H), which separates the A and the W region in Fig. [TJ denotes a Hopf bifurcation 
line. Indeed, within the entire regime left and below of the H line (i.e., not only in the W state), the system exhibits 
oscillatory dynamics corresponding to stable limit cycles (see also [2J, I41H43J ]). Moreover, we find that the Hopf 
bifurcation separating A and W states is supercritical. This implies that, upon crossing the H line, the W state is 
"born" with zero amplitude, but finite frequency at the onset. 

Lowering Ak even further, the character of the oscillatory state (i.e., the limit cycle) changes from wagging to 
tumbling (T). This latter state involves again a periodic motion where the director performs full rotations rather 
than just wagging (W) in the shear plane. We stress, however, that there is no fundamental difference between W 
and T motion in the sense that the change from W to T is not a true bifurcation (see also |23| )• Finally, the curve 
corresponding to the lower boundary of the W/T domain (red line) is a period doubling bifurcation curve (PD1). 
Below the PD1 curve, the W or T state do not correspond to stable limit cycles any more. As seen in the very left of 
Fig. [TJ (i.e., small Ak), there are small deviations between the location of the PD1 curve (obtained by the bifurcation 
analysis) and the colored region of T states obtained by direct numerical integration. These deviations arise from 
the fact that there are very long-living transients in this particular parameter region, making the classification of the 
final state based on the integration scheme rather difficult. However, this "problem" could be fixed by extending the 
numerical integration towards even longer time scales. 



2. Emergence of symmetry-breaking states 

Starting from a W or T state (see upper left area of the diagram Fig. [1]) and reducing the shear rate 7 we cross 
the PD1 curve, thereby entering the domain of kayaking-tumbling (KT) states. Within the KT state, the director 
performs (full) rotations out of the shear plane (contrary to the T state mentioned above). Thus, KT is an example 
of a state breaking the symmetry introduced by the shear flow geometry. 

The nature of the transition separating the W/T and the KT regime depends on the value of the tumbling parameter. 
This is indicated by the solid and dashed parts of the PD1 curve in Fig. [JJ Specifically, in the range Ak ;$ 0.7 (solid 
part), there are two types of stable solutions, that is, the KT solution (a stable limit cycle) and a solution corresponding 
to period-doubled orbits. In practice, however, we never observed these period-doubled orbits, indicating that their 
basin of attraction is much smaller than that of the KT solution. At larger tumbling parameters (Ak <) 0.7), the 
PD1 bifurcation line becomes subcritical. This implies that the period-doubled solutions existing for Ak > 0.7 are 
unstable. Instead, the solution evolves towards the next stable attractor, that is, the KT solution. 

So far we have only considered the PD1 line, below which the T/W limit cycles are non-existent. However, as we 
will discuss in more detail in the subsequent paragraph, the PD1 determines the transition between the T/W and the 
KT regime only when the shear rate is swept from large to small values, that is, when the transition is approached 
from above. If, on the contrary, we increase the shear rate starting from a KT solution, this solution "dies" at a limit 
point of cycle (LPC). These points define the line labeled LPC in Fig. [JJ To complete the discussion of KT states, we 
note that this type of solution can also be reached from the flow-aligned (A) state, provided that the shear rate is not 
too large. This is seen in the right lower corner of the diagram Fig. [TJ The corresponding boundary curve is labeled 
as "LP2" , where LP stands for "limit points" (saddle-node). Crossing this line from the right, KT states are born via 
a "semi-global", that is, a SNIPER or SNIC bifurcation. A SNIC (Saddle-Node-on-an-Invariant Cycle) bifurcation 



results in the appearance of a limit cycle of an infinite period |44j . Indeed, upon approaching the LP2 curve from the 
left (i.e., increasing Ak) the period of the oscillatory motion is found to increase rather steeply (corresponding to a 

-I /Q 

slow-down of the tumbling frequency in a real system). Specifically, the period grows as T ex \b — b c \ , where b is 
the control parameter and b c is the critical parameter value at the SNIPER/SNIC bifurcation J4al46j|. 

Another type of symmetry-breaking state is kayaking- wagging (KW), where the wagging motion of the director 
involves the space outside of the shear plane. As seen from Fig. [JJ the KW state is stable at intermediate values of 
both, the tumbling parameter and the shear rate. One way to reach the KW regime is to start from the W region 
(at suitable values of 7, i.e. 2.4 < 7 < 4.8) and to increase Ak- Doing this, one crosses a line of branch points of 
cycles (BPC). Beyond this line, the (in-plane) wagging (W) ceases to exist, and a limit cycle of type KW is born. We 
note that the appearance of BPCs implies that the KW oscillations are born in a pairwise manner. Thus, there is 
always a coexistence of two KW solutions in the parameter region labeled as KW in Fig. [TJ Upon increasing Ak even 
further, the KW disappears. As revealed from Fig. [TJ this happens either at a second, subcritical period doubling 
curve (labeled PD2) or via a Ncimarck-Sackcr (NS) bifurcation. The latter leads into a regime of chaotic solutions 
which are not further considered here (for a discussion of such states, see |28j). 



3. Bistability 

As already mentioned in the preceding paragraph, there are parameter regions where the dynamical behavior of 
the sheared system is characterized by bistability (we do not further consider here the trivial coexistence of two KW 
states within the KW region) To find bistable regions, we have investigated whether and how the results of the (direct) 
numerical integration change when the parameters are varied in a different way. The resulting bistable regions are 
visualized in the two parts of Fig.[2j The first region is bounded by the LPC line (from above), the dashed (subcritical) 
part of the PD1 line (from below), and the BPC line (from the right). Depending on whether one enters that area 
from high or low shear rates, either in-plane T/W oscillations or out-of-plane KT oscillations are found. The second 
region is again bounded by the LPC line from above, whereas the left and right boundaries are the line of BPCs and 
the curve labeled PD2, respectively. Given the existence of these bistable regions it would be very interesting to see 
the consequences of this "coexistence" of dynamical states, if we allowed the system to be inhomogeneous. Indeed, 
first investigations in this direction [47] already revealed the appearance of spatial domains characterized by different 
orientational states. These domains interact and change with time; an example being the growth of one domain at 
the expense of another one. A detailed study of an inhomogeneous system is under way, but is outside the scope of 
the present paper. We note in this context that there are already some studies investigating inhomogeneous systems. 
In particular, Das et. al [48( studies, based on the same model we are using, inhomogeneous systems for similar shear 
rates and coupling parameters. However, the main focus was on detecting spatio-temporal chaos. 

B. Beyond the nematic phase 

So far we have discussed the shear-induced nonlinear dynamics at the temperature 9 — corresponding to the 
nematic phase of the equilibrium system. Given the variety of dynamic states observed at this low temperature, it is 
tempting to investigate the "fate" of these dynamical states when 9 is increased towards values where the equilibrium 
state becomes isotropic. Clearly, under isotropic equilibrium conditions, any orientational ordering is induced by the 
shear flow. Indeed, previous theoretical investigations have already shown that shear (applied to an originally isotropic 
system) can generate a flow- aligned state, paranematic state; moreover, it can shift the transition into the nematic 
state jg. 

However, usually this shear-induced, paranematic ordering is assumed to be stationary. This is also seen exper- 
imentally, e.g., for fd-viruses at concentrations below the paranematic- nematic transition [18] . We note, however, 
that the particles involved in these experiments correspond to particular values of the tumbling parameter, Ak- For 
example, as described below Eq. (|2.11[) . fd-viruses are characterized by tumbling parameters of about Ak ~ 0.7 — 0.8. 
There remains the question whether oscillatory dynamic states could exist beyond the nematic phase, if the particle's 
properties and thus, the tumbling parameter Ak is changed. 

To explore this issue within our mesoscopic approach we have extended the calculations at 9 = described in 
Sec. IIII Al towards a range of larger values of 9, corresponding to higher real temperatures (smaller real concentrations) 
in thermotropic (lyotropic) systems. An exemplary state diagram in the Ak-7 plane is shown in Fig. [3J The reduced 
temperature is fixed to 9 = 1.20, which is well above the value 9 = 9/8 = 1.125 (at T u or c u , respectively) corresponding 
to the upper limit of the (metastable) nematic phase. The curves and marked points have been obtained by MATCONT 
[39i |. In addition to these bifurcation lines, we have indicated in Fig.[3jthc degree of nematic order induced by the shear 
flow. Specifically, the colors are chosen according to the value of the eigenvalue ^ that characterizes the ordering in 
the direction of the director. In other words, the eigenvalue /i^ characterizes the degree of uniaxial ordering. Note, 
however, that at nearly all parameter combinations shown in Fig. [3] the orientational order is actually biaxial. 

From Fig. [3J we see that, at small shear rates and tumbling parameters, the system is only weakly ordered, as 
expected at the high temperature (low concentration) considered. Increase of either Ak or 7 then leads to an increase 
of the order and hence larger /Z3, reflecting shear-induced alignment. 

However, the most important information from Fig.[3]is that there is a range of shear rates and tumbling parameters 
where the shear induces oscillatory motion. This parameter range is represented by the "island" in the left part of 
the figure, corresponding to tumbling parameters in the range Ak ^ 0.75. Within this island, the director performs a 
wagging (W) motion within the shear plane, implying that the degree of alignment oscillates in time. This is indicated 
by the grainy color inside the W region. The grainy colors stem from different values of the largest eigenvalue /X3 
during the oscillation. For every pair of parameters we determine 113 only in one point of time. The upper boundary 
of the oscillatory island is a Hopf bifurcation curve, whereas the lower boundary is a homoclinic orbit. 
To illustrate the behavior of the largest eigenvalue /13 of the alignment tensor as the order parameter versus the shear 
rate in this parameter region we provide a codimension-1 bifurcation analysis in Fig.2]for fixed Xk = 0.7 (see thin red 
line in Fig. [3]). At the limit points (LP) and at the Hopf point (H) the stability of the fixpoint changes. The dashed 
lines indicate an unstable fixpoint. The red lines depicts the minima and maxima of the limit cycle (in the projection 





FIG. 2. (Color online) Enlarged views of the parameter section of Fig. [T] in which bistability occurs. The data in part (a) were 
obtained by sweeping the shear rate from higher to lower values, whereas in part (b), the shear rate was swept from lower to 
higher values. In each step, the final values of the integration is used as an initial condition for the next integration step. To 
avoid staying on unstable solutions, the initial condition was perturbed by a small noise. 
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FIG. 3. (Color online) State diagram at the reduced temperature 9 = 1.20, where the equilibrium state (7 = 0) is isotropic. 
The curves are obtained from a codimension-2-bifurcation analysis. The colors here reflect the value of the eigenvalue fi = (13, 
which may be considered as the degree of the uniaxial ordering (see color bar on the right side of the figure) . One observes an 
"island" of wagging (W) states. Region Al corresponds to stable, shear- induced nematic order (flow-alignment). In region A2, 
there is a coexistence of a paranematic and a nematic state. A closer view of this region (revealing even more states) is given 
in Fig. [5] The thin red line indicates the path of the codimension-1 bifurcation, for fixed Xk = 0.7, that is discussed in Fig. [4] 



to the largest eigenvalue 113) that is evolving from the Hopf bifurcation point (H). The period T of the limit cycles 
tends to infinity, when the parameter is close to the homoclinic bifurcation (Horn) (colored brown). The period grows 
logarithmically as the bifurcation is approached, i.e. T ex l og( \b — b c \), where b is the control parameter and b c is the 
critical parameter value at the homoclinic bifurcation [45|, |4y, |49( ■ Also note the hysteresis taking place in the small 
parameter region between the lower saddle-node bifurcation (LP) and the homoclinic bifurcation (Horn). 
To our knowledge, this is the first time that oscillatory states outside of the nematic region of the (equilibrium) 
phase diagram are found. One reason might be that most earlier theoretical studies based on the present, mesoscopic 
approach [38l . |48[ focus on larger values of the tumbling parameter. Indeed, in the range Ak <] 0.75 the present 
calculations reproduce the (well-established) fact of a shear-induced, first-order, paranematic-nematic transition [9U38J. 
Consider, as an example, the case Ak = 0.8. Increasing the shear rate from small values, the system remains essentially 
disordered (as indicated by the blue color) up to 7 ps 0.2. At this point, one crosses the LP bifurcation line and enters 
a bistable region characterized by two stable fixed points, that correspond to a paranematic ordering coexisting with 
a (stationary) shear-aligned state. In a plot of the order parameter versus shear rate (not shown here), this bistable 
region (labeled A2 in Fig. [3]) would correspond to the hysteresis region. Finally, upon increasing 7 even further, the 
pscudo nematic ordering becomes unstable and one enters the region Al. Here, the order parameter is relatively large, 
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FIG. 4. (Color online) Codimension-1 bifurcation diagram with the shear rate j as control parameter and the largest eigenvalue 
/i = H3 of the alignment tensor as the order parameter. The analysis is performed along the thin red line shown in Fig. [3] 
and Fig. [5] (9 = 1.20, Xk = 0.7, a = 0). The solid lines depict stable fixpoints. Dashed lines indicate unstable fixpoints. 
Between the homoclinic bifurcation (Horn, colored brown) and the supercritical Hopf bifurcation (H), there are stable limit 
cycles corresponding to the wagging state (W). These limit cycles are illustrated by vertical lines between the maximal and 
minimal values of the limit cycles for the projection on the largest eigenvalue fi = /13. 



reflecting pronounced shear-alignment. 

The bifurcation scenario occurring at the parameters where the wagging island, on the one hand, and the A2/A1 
region, on the other hand, approach one another, is rather complicated. This is illustrated in Fig. [5J which shows an 
enlarged view of the relevant section of Fig. [3] 

In particular, one can see that there is a small, additional region labeled A3. In this region one finds both, a 
stable fixed point representing weak flow- alignment, and a stable limit cycle representing dynamical wagging. In 
other words, A3 is again a bistable region. The upper boundary of this region is either a limit point (LP) saddle-node 
bifurcation curve (smaller values of Ak) or a Hopf bifurcation curve (larger values of Ak)- Finally, the lower boundary 
is a homoclinic curve. Below the homoclinic curve the wagging solutions does not exist any more. 
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FIG. 5. (Color online) Enlarged view of the parameter section of Fig. [3] at the border of the wagging island and the A2/A1 
regions (6 — 1.20). Within the region labeled A3, stationary alignment (with low value of the order parameter) coexists with 
an oscillatory (W) solution. 



C. Non-equilibrium phase diagrams 



In the preceding paragraph we have studied the system in the parameter plane spanned by the shear rate, 7, and 
the coupling parameter Ak, the latter being related to the shape of the particles [see Eq. (|2.11|) ]. Thus, it is clear 
that Ak is not an adjustable quantity in a specific, experimental system. Here, the relevant parameters are rather the 
temperature or concentration, respectively, and the shear rate or its conjugate, the shear stress, while the tumbling 
parameter Ak is constant (up to polydispcrsity effects). 

Motivated by this experimental situation we present in Figs. [6][8] state diagrams in the 8-j or 07 plane. These 
diagrams may be considered as non-equilibrium phase diagrams (as a generalization of the equilibrium phase diagrams 
corresponding to the case 7 = 0). We have obtained these diagrams for three values of the tumbling parameter Ak, 
that is, at Ak = 1.25, 0.74, and 0.6. These values seem particularly interesting judging from our previous analysis in 

Sees. ImATnTBl 

Iii the subsequent figures [B][51 the lower x-axis shows the concentration, c, while the upper 2- axis shows additionally 
the reduced temperature. The y-axis denotes the shear rate. By choosing c as a main variable, we take into account 
the fact that many recent experimental studies [17|, [HI focus on lyotropic systems. To map the quantities c and 9 
onto each other we have employed the second member of Eq. (|2.3[) that can be solved with respect to c. Clearly, 
application of the resulting equation requires an input for the clearing point concentration, ck, and the pseudo-critical 
concentration, c*. Here we have arbitrarily chosen ck = 0.4 and c* = 0.5. 

In Fig. [51 we have indicated ck (corresponding to 6 = 1) and c* (8 = 0) by colored vertical lines. In addition, we 



13 




FIG. 6. (Color online) Non-equilibrium phase diagram in the plane spanned by the concentration (bottom horizontal axis) 
or temperature (top horizontal axis) and the shear rate (vertical axis) at fixed tumbling parameter Ak = 1.25. The colors 
represent the eigenvalue fi — /13 of the alignment tensor a. We also show the bifurcation lines obtained by a codimension-2 
analysis. In the triangular area labeled A2, we find bistability between a paranematic and a nematic state with flow-alignment. 



have put a vertical line at the concentration c u — 0.3875 (9 = 9/8 = 1.125) below (above) which the nematic state 
is absolutely unstable. Finally, the color of the different regions in Fig. [6] indicate the degree of nematic order [see 
Eq. (|2.6|) ]. neglecting biaxiality. Thus, the lower x-axis corresponds to the equilibrium phase diagram of a lyotropic 
system. 

We start by considering the case Ak = 1-25 illustrated in Fig. [51 At small concentrations below the stability limit of 
the nematic phase (i.e., left from the red vertical line) and small, yet non-zero, values of 7 we observe a paranematic 
state with very weak nematic order. Increasing the shear rate at a fixed c in this range of concentrations, the system 
enters first a region (labeled A2) of bistability between the paranematic and a shear-aligned state, and finally the true 
shear-aligned state (Al). These phenomena are characteristic of the well-known, shear-induced paranematic-nematic 
transition already mentioned in Sec. lIIIBl Indeed, following earlier studies we can interpret the two violet lines labeled 
by LP1 as "spinodals" of the isotropic-nematic transition under shear and their merging point (c c , 7 C ) as "critical 
point" of this transition. The analogy to an equilibrium critical point (of, e.g. the vapor-liquid transition of a fluid) is 
that at shear rates above the critical point, the transition becomes continuous, that is, the order parameter increases 
smoothly when c is increases at fixed 7 > -%. Consistent with earlier studies [9|, [5(|, we find that the spinodal lines 
are curved towards the left. This is plausible, since shear tends to enhance the degree of nematic ordering and thus 
supports the transition into a (shear-) aligned state. We note that similar effects in the non-equilibrium phase diagram 
have been observed in studies of the impact of biaxial stretching flow [25[ . 

We also find from Fig. [5] that at Ak = 1.25, oscillatory states only appear at high concentrations (low 9) within the 
nematic regime. Consistent with our analysis of the case 9 = 0/c = c* in Sec. MI Al the states are of kayaking-tumbling 
(KT) type in the range of shear rates considered (compare to Fig.QJ. 

In Sec. IIII Bl we have found that at smaller values of Ak, oscillatory states can become stable even at reduced 
temperatures or concentrations where the equilibrium system is isotropic. Motivated by this observation, we present 
in Figs. [7] and |S] 07 diagrams at Ak =0.74 and Ak = 0.6. 

Similar to the case of Ak = 1.25, one may identify the spinodal of the paranematic-nematic transition, with the 
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FIG. 7. (Color online) Same as Fig. [6] but for Ak = 0.74. In the area right of the Hopf bifurcation curve (H) we observe 
oscillatory states. Specifically, these are Wagging/Tumbling (W/T) in the area between the Hopf (H) and the period doubling 
(PD) bifurcation curve, and Kayaking/Tumbling (KT) in the area right of the PD curve. In between the LPC and the PD 
curve, we find bistability between these states. 



"critical" shear rate 7 C being somewhat larger than before. This mean that higher shear rates are necessary to see 
a complete disappearance of the first-order transition. However, the most important new feature appearing at lower 
values of Ak is that the Hopf (H) curve, which marks the onset of oscillatory states, moves (at least partially) into 
the parameter region left of the red, vertical line, that is, into the regime where the equilibrium ncmatic state is 
absolutely unstable. These oscillatory states are of type wagging (W) or tumbling (T) [as argued before, there is 
no substantial difference between these motions from the perspective of the bifurcation analysis]. Increasing the 
concentration further, systems at Ak = 0.74 or Ak = 0.6 first enter a region where, depending on initial conditions, 
the dynamics is either of type W/T or of type KT. This region is located in between the LPC curve (limit point of 
cycle) and the PD (period doubling) curve. Finally, at concentrations beyond the PD curve, only T states are stable. 

Coming back to the regime of low concentrations we see from Fig. [7] that at Ak = 0.74, the oscillatory region left 
of the (nematic) stability limit is rather small; in fact, the major part of the corresponding Hopf line lies inside the 
spinodal. Interpreted in terms of a real system, this would mean that oscillations do occur, but only at concentrations 
larger than the critical one (defined by the top of the spinodal), i.e., inside the true nematic phase of the sheared system. 
This is indeed consistent with recent experimental results for fd- viruses under shear flow. In these experiments, the full 
paranematic- nematic binodal at shear rates 7 > was mapped out. The "coexisting" nematic phase is characterized 
by tumbling motion, whereas the paranematic state is flow-aligned. The tumbling-to-aligning transition line then 
ends at the maximum of the binodal (la , [la . l5l| . 

The behavior just discussed changes at even smaller values of the tumbling parameter such as the case Ak = 0.6. 
Indeed it is seen from Fig. [5] that, upon increasing 7 from zero in the relevant parameter region (i.e., left of the red, 
vertical line), the Hopf curve starts only at the critical point (7 c ,c c ) and then initially bends to the left. In other 
words, W or T states can exist already at the less ordered side of the paranematic-nematic transition, consistent with 
the 7-Ak diagram discussed in Sec. IIIIBI [sec Fig. [3]. The corresponding behavior of the order parameter /X3, i.e., 
the largest eigenvalue of the alignment tensor, as function of c (9) is shown in Fig. O Only at larger values of the 
shear rate we finally observe a bending of the Hopf curve towards higher concentrations. This latter behavior may be 
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FIG. 8. (Color online) Non-equilibrium phase diagram for Ak = 0.6. Compared to Fig. [7](Ak = 0.74), the wagging region (W) 
is more extended towards the concentration (temperature) region, where the equilibrium state is isotropic (i.e. left from the 
red vertical line). The thin white horizontal line indicates the path for the codimension-1 diagram presented in Fig. [9] 

interpreted such that strong shear typically favors shear- aligned states (Al) rather than oscillatory motion (compare 
with Fig.l). 



IV. CONCLUSIONS 



In the present study we have employed a numerical path continuation analysis combined with mesoscopic equations 
of motion to investigate the complex dynamical behavior of a homogeneous system of rod-like particles under planar 
(Couette) shear flow. Compared to a conventional integration of the mesoscopic equations [l{| [20, [29|, bMJ , a main 
advantage of the continuation method is that it provides not only full information about the (long-time) dynamics at 
different system parameters, but also about the nature of the bifurcation lines separating different regions of the non- 
equilibrium phase diagram. This information is not only of theoretical, but also of practical interest. For example, 
close to the bifurcation line separating different states, the continuation method can predict the behavior of the 
amplitude and frequency of an order parameters when one crosses the bifurcation line by varying a suitable system 
parameter (e.g., shear rate, concentration, . . . ). Furthermore, the method proved to be very helpful in determining 
multistable regions, that might otherwise be missed due to their small extension in parameter space. On top of these 
insights, however, we have found a new aspect of the non-equilibrium behavior which was formerly unknown. Namely, 
the sheared system can exhibit an oscillatory state (wagging) at temperatures/concentrations outside the nematic 
region of the equilibrium phase diagram. To our knowledge this is the first time not only in theory, but also from 
the experimental perspective that oscillatory states outside the nematic region have been observed. One possible 
reason might be that, according to our results, the oscillatory states only occur for rather small values of the coupling 
parameter Xk, corresponding to rod- like particles with relatively small aspect ratio. Real systems such as suspensions 
of fd- viruses typically involve more elongated particles, for which our theory predicts oscillatory states only within the 
nematic region, consistent with experiments [151 ] . Indeed, our non-equilibrium phase diagram (obtained by performing 
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FIG. 9. (Color online) Codimension-1 diagram using the bifurcation parameter c (9). The path of the bifurcation analysis is 
along the white horizontal line in Fig. [8] (7 = 0.9, Ak = 0.6, a = 0). Solid (dashed) lines indicate stable (unstable) fixpoints. 
The limit cycles are illustrated by vertical lines between the maximal and minimal values for the projection on the largest 
eigenvalue \i — /i3 as the order parameter. Note the bistable region of wagging states (red limit cycles) and kayaking tumbling 
states (blue limit cycles) between the period doubling bifurcation (PD) and the limit point of cycle (LPC) bifurcation. 



51|. 



calculations for a large range of concentrations) has strong similarities with that presented in Rcf. 

Based on the present results one can now proceed towards the investigation of inhomogencous systems, including the 
question how different spatial regions interact within the multistable regions of the parameter space. A particularly 
interesting aspect in this context concerns the front speeds [52, [53[ and the basins of attraction of the involved states. 
Another motivation of including inhomogeneities is to explore the appearance of "banded" states, such as vorticity 
banding (which has already been observed in experiments |54J ) and gradient banding [5l| . Work in these directions 
is in progress. 



Appendix A 



In this Appendix we present some considerations regarding the structure and the symmetries of Eqs. 
To start with, we rewrite the right-hand sides of these equations. The Landau-de Gennes free energy in Eq. (|2.4[) can 



be written in projections of the tensorial basis [37[ as 



$ = \{e + o 2 )o 2 - a 
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2 2 -3(al + a 2 2 ) + -(a 2 3 + al) 
-W?> [ai{a\ - al) + 2a 2 a s a 4 \ . (Al) 

InEqs. flAT]), a 2 = ^ =Q a 2 . 

We now consider several cases. 

i) For vanishing shear rate 7 = the dynamic equations reduce to 

4-^ = -$, = -^-$- (A2) 

at oat 

The $i are identical to those appearing in Eqs. (|2.13[) . From Eq. (|A2[) we see that the alignment follows a gradient 
dynamics. Therefore, no oscillatory solutions are possible for 7 = 0, i.e., oscillatory solutions are all flow-induced |55j . 
To see this, we determine the Jacobian J of the right hand side of Eq. (|A2[) as 

d d d d d 

Jij = ^~(-*i) = —^-(tt®) = -^-(^) = J n- ( A3 ) 

aaj aaj aai aai aaj 

The Jacobian «Zy is symmetric due to the symmetry of the second-order derivatives. A symmetric matrix has only 
real eigenvalues. Thus, when performing a linear stability analysis based on «/, we cannot find any oscillations, 
ii) For the case 7 7^ but a = 0, it is instructive to introduce a function H defined as 

H = ^(a 2 + a 2 2 + ^(at + al)-V3\ K a 1 ). (A4) 

Note that H is independent of ao and proportional to the shear rate 7. With this function, Eqs. (|2. 12112. 13]) can be 
rewritten as 



d 
dT = 


_A$ 

da 


d 
dt ai = 


d d 

dai da 2 


d 
dT 2 = 


oa 2 aa\ 


d 
dt a3 = 


9 $ 9 if 

9fl3 9a4 


d 
dT 4 = 


9a4 day, 



(A5) 

The structure of the Eqs. (|A5[) reflects the presence of two oscillatory subunits (01,02) and (03,04). The partial 
derivative of H with respect to a 2 appears in the dynamical equation for ai and vice versa the partial derivative of H 
with respect to ai in the dynamical equation for a 2 . Similarly, the partial derivative of H with respect to 04 appears in 
the dynamical equation for 03 and vice versa the partial derivative of H with respect to 03 in the dynamical equation 
for 04. In particular, oscillations in a\,a 2 with 03 = 04 = indicate tumbling or wagging within the shear plane. If 
there are additionally oscillations in 03 and 04 we have out of shear plane kayaking states. 

iii) In order to make a similar reformulation of the dynamics for the more general case 7 7^ and a 7^ 0, two additional 
dissipative couplings between ao and a 2 as well as between 03 and 04 must be added. In this spirit, we amend the 
potential <J? in Eqs. (|A5[) by an additional term V (i.e., $ — > $ + V). For the amendment we have 

v = a i(y a 3ai - -^V3a a 2 ). (A6) 

We note that the additional term V for a 7^ only contributes to the relaxational part of the dynamic equation. 
So in general the Eqs. (|2.12j|2.1^|) can be expressed as derivatives of $, V and H . 
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Appendix B 

In this Appendix we give a short description on how we determined our bifurcation diagrams for homogeneous 
systems based on the freely available software package MATCONT [39|. As shown in Sec. Ill Bt the present dynamical 
system can be formulated as dx/dt = f(x, (3), where the vector x contains the components a^,ai,. . . ,04 and (3 contains 
the bifurcation parameters. Here we are dealing with a codimension-2 diagram, where one has two bifurcation 
parameters (all other parameters are fixed). We investigate mainly two cases. In the first case the bifurcation 
parameters are Ak and 7, whereas the dimensionless temperature 9 [see Eq. (|2.3j) ] is kept fixed. 

If we now choose a set of parameters (3 = (3q and start to integrate the system with initial conditions x = Xo, there 
are three possibilities for the outcome of this procedure (disregarding transients). First, the solution goes to a stable 
fixed point (equilibrium point/EP); second, it goes to a stable periodic oscillatory solution (limit cycle/LC) or third, 
it goes to a possibly irregular and chaotic attractor. From a practical point of view, it is most convenient to choose 
Xo and /3o such that one ends up in an EP or LC solution. Starting from that solution one changes one parameter, 
i.e., one performs a codimcnsion-1 bifurcation analysis. 

The continuer algorithm then tries to stay on the EP or LC solution, while also monitoring the Jacobian of 
the system and its eigenvalues. Typically one can encounter saddle node or limit point bifurcation (LP), a Hopf 
bifurcation H, a limit point of cycles (LPC), a torus or Neimarck-Sacker bifurcation (NS), and a flip or period 
doubling bifurcation (PD) (see |5g). These codimension-1 bifurcations points are then used as starting points for a 
codimension-2 bifurcation analysis. In this way one can calculate branches of LP, H, LPC, NS or PD depending on 
two bifurcation parameters. Following these branches one can also encounter special codimension-2 points like a Cusp 
point (CP), Bogdanov-Takens point (BT), Zero Hopf (ZH) or Generalized Hopf (GH). We note that this list is by no 
means complete. 

A characteristic feature of codimension-2 points is that they correspond to an intersection of different branches of 
bifurcations. For example, a BT corresponds to the intersection of a Hopf branch, a LP branch and a homoclinic 
branch. Based on that knowledge one can now attempt to continue along the missing branches. Reiteration of these 
steps finally yields the full codimcnsion-2 bifurcation diagram. 

We are noting that, within our bifurcation analysis, both KT and KW solutions are indeed rather hard to find. The 
reason is that these symmetry-breaking states are not born in local bifurcations (such as the Hopf bifurcation separat- 
ing A and W states) but rather in global bifurcations. To identify such a bifurcation it is not sufficient to just monitor 
the Jacobian and its eigenvalues. However, once one has obtained a limit cycle solution of type KW or KT, the nu- 
merical continuation via MATCONT allows to follow these solutions in parameter space, until the boundary is reached. 
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